# Figure 4: lecture 24
set.seed(10)
#generate data
stuff<-c(rnorm(1000), rnorm(300,mean=5.5))
# kernel density estimate
out.d <- density(stuff)
# Determine W that cuts off p = 0.95
out.diff <- sapply(seq(0.35, 0.5,.001), function(x) sum(out.d$y[out.d$y>=quantile(out.d$y, x)])/sum(out.d$y)-0.95)
which.one <- which.min(out.diff[out.diff>0])
cut.off <- seq(0.35,0.5,.001)[which.one]

# credible intervals
my.x<-out.d$x[out.d$y>=quantile(out.d$y, cut.off)]
my.y<-out.d$y[out.d$y>=quantile(out.d$y, cut.off)]
my.x2 <- out.d$x[out.d$x<=quantile(stuff,.975) & out.d$x>=quantile(stuff,.025)]

# identify regions under curve
grp1.x <- my.x[my.x >=range(my.x[my.x<3])[1] & my.x <=range(my.x[my.x<3])[2]]
grp1.y <- my.y[my.x >=range(my.x[my.x<3])[1] & my.x <=range(my.x[my.x<3])[2]]
grp2.x <- my.x[my.x >=range(my.x[my.x>3])[1] & my.x <=range(my.x[my.x>3])[2]]
grp2.y <- my.y[my.x >=range(my.x[my.x>3])[1] & my.x <=range(my.x[my.x>3])[2]]

plot(density(stuff), xlab=expression(theta), main='',cex.lab=1.2, col='grey70')
rug(my.x,side=1,col='dodgerblue1',line=-.5)
rug(my.x2,side=1,col=2)
#shade area under curve
polygon(c(grp1.x, grp1.x[length(grp1.x)],grp1.x[1]),c(grp1.y,0,0), col='lightblue', border=NA, density=50)
polygon(c(grp2.x, grp2.x[length(grp2.x)],grp2.x[1]),c(grp2.y,0,0), col='lightblue', border=NA, density=50)
segments(grp1.x[1],quantile(out.d$y, cut.off),grp2.x[length(grp2.x)],quantile(out.d$y, cut.off),lty=2)

arrows(3,.15,0.5,.1,code=2,angle=45,length=.1)
arrows(3.3,.15,5.5,.07,code=2,angle=45,length=.1)
text( 3.15,.15, 'Area = 0.95',pos=3,offset=0.1)
legend('topright',c('HPD','Percentile'),col=c('dodgerblue1',2), pch=15, pt.cex=1.2, cex=.9, bty='n', title=expression(bold('95% Credible Intervals')))
lines(out.d$x, out.d$y)
text(grp2.x[length(grp2.x)],quantile(out.d$y, cut.off),'W',pos=4)